Sister pair DE analysis, doplots and volcanoplots ctrl -vs- lumb
Fabio Sacher
18.06.2024
library(Seurat)## Attaching SeuratObjectlibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(stringr)
library(tibble)
library(patchwork)
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutDE table
First we load the sister pair DE tables and filter for:
absolute avg_log2FC > 0.5 (~41% increase)
p_val_adj < 0.01
DE_list <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_lumb_sis_markers.rds")
for (i in seq(DE_list)) {
DE_list[[i]] <- DE_list[[i]] %>%
arrange(desc(avg_log2FC)) %>%
filter(abs(avg_log2FC) > 0.5) %>%
filter(p_val_adj < 0.01)
}
DE_table <- do.call(rbind, DE_list)
dim(DE_table)## [1] 985 8delta pct distribution
par(mfrow = c(2,2))
hist(abs(DE_list[[1]]$delta_pct), breaks = 20)
abline(v = 0.1, lty = "dashed", col = "red")
hist(abs(DE_list[[2]]$delta_pct), breaks = 20)
abline(v = 0.1, lty = "dashed", col = "red")
hist(abs(DE_list[[4]]$delta_pct), breaks = 20)
abline(v = 0.1, lty = "dashed", col = "red")
hist(abs(DE_list[[5]]$delta_pct), breaks = 20)
abline(v = 0.1, lty = "dashed", col = "red")
Now we filter the DE lists for absolute delta percentage > 0.1.
for (i in seq(DE_list)) {
DE_list[[i]] <- DE_list[[i]] %>%
filter(abs(delta_pct) > 0.1)
}
DE_table <- do.call(rbind, DE_list)
dim(DE_table)## [1] 778 8Broad clusters
broad_order <- c("progenitors",
"FP",
"RP",
"FP/RP",
"neurons",
"OPC",
"MFOL",
"pericytes",
"microglia",
"blood",
"vasculature"
)Integrated data
Load the integrated control and poly data.
int_path <- "Gg_ctrl_lumb_int_seurat_250723"
my.se <- readRDS(paste0("~/spinal_cord_paper/data/", int_path, ".rds"))
annot_int <- read.csv(list.files("~/spinal_cord_paper/annotations",
pattern = str_remove(int_path, "_seurat_\\d{6}"),
full.names = TRUE))
if(length(table(annot_int$number)) != length(table(my.se$seurat_clusters))) {
stop("Number of clusters must be identical!")
}
# rename for left join
annot_int <- annot_int %>%
mutate(fine = paste(fine, number, sep = "_")) %>%
mutate(number = factor(number, levels = 1:nrow(annot_int))) %>%
rename(seurat_clusters = number)
ord_levels <- annot_int$fine[order(match(annot_int$broad, broad_order))]
# add cluster annotation to meta data
my.se@meta.data <- my.se@meta.data %>%
rownames_to_column("rowname") %>%
left_join(annot_int, by = "seurat_clusters") %>%
mutate(fine = factor(fine, levels = ord_levels)) %>%
mutate(seurat_clusters = factor(seurat_clusters, levels = str_extract(ord_levels, "\\d{1,2}$"))) %>%
column_to_rownames("rowname")
ctrl_poly_int_combined_labels <- readRDS("~/spinal_cord_paper/annotations/ctrl_lumb_int_combined_labels.rds")
my.se <- AddMetaData(my.se, ctrl_poly_int_combined_labels)DimPlot
DimPlot(
my.se,
group.by = "annot_sample",
reduction = "tsne",
label = TRUE,
repel = TRUE
) +
NoLegend()## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Cluster order
Get the cluster order from the spearman correlation heatmap of the control and poly integrated data. Then we filter for the neuronal clusters only.
corr_heatmap <- readRDS("~/spinal_cord_paper/output/heatmap_spearman_ctrl_lumb.rds")
#heatmap order
htmp_order <- data.frame("label" = corr_heatmap[["gtable"]]$grobs[[4]]$label) %>%
mutate(label = str_remove(label, "_int")) %>%
mutate(label_ordered = paste(str_sub(label,6 ,-1), str_sub(label, 1, 4), sep = "_"))
my.se@meta.data <- my.se@meta.data %>%
mutate(annot_sample = factor(annot_sample, levels = htmp_order$label_ordered))
Idents(my.se) <- "annot_sample"
# filter for the neuronal clusters
my.se <- subset(my.se, idents = htmp_order$label_ordered[grepl("neurons|MN", htmp_order$label_ordered)])
DimPlot(
my.se,
group.by = "annot_sample",
reduction = "tsne",
label = TRUE,
repel = TRUE
) +
NoLegend()
my.se@active.assay <- "RNA"Individual dot plots
# select top50 by log2FC
for (i in seq(DE_list)) {
DE_list[[i]] <- DE_list[[i]] %>%
slice_max(order_by = abs(avg_log2FC), n = 50) %>%
arrange(desc(avg_log2FC))
}
p1 <- modplots::mDotPlot2(my.se,
group.by = "annot_sample",
assay = "RNA",
# reverse order of DE genes so number one is on top
features = rev(DE_list[[1]]$Gene.stable.ID),
gnames = modplots::gnames,
cols = c("lightgrey", "black")) +
theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
coord_flip() +
xlab(names(DE_list)[1])## p2 <- modplots::mDotPlot2(my.se,
group.by = "annot_sample",
assay = "RNA",
# reverse order of DE genes so number one is on top
features = rev(DE_list[[2]]$Gene.stable.ID),
gnames = modplots::gnames,
cols = c("lightgrey", "black")) +
theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
coord_flip() +
xlab(names(DE_list)[2])
p3 <- modplots::mDotPlot2(my.se,
group.by = "annot_sample",
assay = "RNA",
# reverse order of DE genes so number one is on top
features = rev(DE_list[[3]]$Gene.stable.ID),
gnames = modplots::gnames,
cols = c("lightgrey", "black")) +
theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
coord_flip() +
xlab(names(DE_list)[3])
p4 <- modplots::mDotPlot2(my.se,
group.by = "annot_sample",
assay = "RNA",
# reverse order of DE genes so number one is on top
features = rev(DE_list[[4]]$Gene.stable.ID),
gnames = modplots::gnames,
cols = c("lightgrey", "black")) +
theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
coord_flip() +
xlab(names(DE_list)[4])
p5 <- modplots::mDotPlot2(my.se,
group.by = "annot_sample",
assay = "RNA",
# reverse order of DE genes so number one is on top
features = rev(DE_list[[5]]$Gene.stable.ID),
gnames = modplots::gnames,
cols = c("lightgrey", "black")) +
theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
coord_flip() +
xlab(names(DE_list)[5])pdf("~/spinal_cord_paper/figures/ctrl_lumb_dotplot_individual.pdf", height = 13, width = 20)
(p1 + p2 + p3 + p4 + p5) + plot_layout(guides = "collect", nrow = 1)
dev.off()## png
## 2Volcanoplots
p.adj <- 0.01
l2fc <- 0
# select top50 by log2FC
for (i in seq(DE_list)) {
DE_list[[i]] <- DE_list[[i]] %>%
mutate(delta_pct_sign = case_when(
delta_pct < 0 ~ "-",
delta_pct > 0 ~ "+",
delta_pct == 0 ~ "0"
))
}
toplot <- do.call(rbind, DE_list) %>%
rownames_to_column("contrast") %>%
mutate(contrast = str_remove(contrast, "\\.\\d{1,2}")) %>%
mutate(contrast = str_replace_all(contrast, " ", "_"))
volplot <- ggplot(data = toplot,
aes(x = avg_log2FC,
y = -log10(p_val_adj),
label = Gene.name,
color = delta_pct_sign,
size = abs(delta_pct)
)) +
geom_point(shape = 21) +
geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
scale_color_manual(values = c("#419c73", "black")) +
scale_size_continuous(range = c(0.5, 4)) +
facet_wrap("contrast", ncol = 5) +
ylab("-log10(padj)") +
theme_bw()
ggplotly(volplot)pdf("~/spinal_cord_paper/figures/Fig_4_volcanoplots.pdf", width = 15, height = 5)
(volplot +
ggrepel::geom_text_repel(size = 3, color = "black"))## Warning: ggrepel: 43 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: ggrepel: 44 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: ggrepel: 42 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: ggrepel: 13 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps## Warning: ggrepel: 21 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps# Date and time of Rendering
Sys.time()## [1] "2024-06-18 10:23:50 CEST"sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS/LAPACK: /scicore/soft/apps/OpenBLAS/0.3.1-GCC-7.3.0-2.30/lib/libopenblas_sandybridgep-r0.3.1.so
##
## locale:
## [1] en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] plotly_4.10.0 patchwork_1.1.1 tibble_3.1.8 stringr_1.4.0
## [5] ggplot2_3.3.3 dplyr_1.0.10 SeuratObject_4.0.2 Seurat_4.0.5
##
## loaded via a namespace (and not attached):
## [1] plyr_1.8.6 igraph_1.2.6
## [3] lazyeval_0.2.2 sp_1.4-5
## [5] splines_4.1.0 crosstalk_1.1.1
## [7] listenv_0.8.0 scattermore_0.7
## [9] GenomeInfoDb_1.28.0 digest_0.6.27
## [11] htmltools_0.5.1.1 fansi_0.5.0
## [13] magrittr_2.0.1 memoise_2.0.0
## [15] tensor_1.5 cluster_2.1.2
## [17] ROCR_1.0-11 globals_0.16.2
## [19] Biostrings_2.60.0 matrixStats_0.58.0
## [21] modplots_1.0.0 spatstat.sparse_3.0-0
## [23] colorspace_2.0-1 blob_1.2.1
## [25] ggrepel_0.9.1 xfun_0.34
## [27] RCurl_1.98-1.3 crayon_1.4.1
## [29] jsonlite_1.7.2 spatstat.data_3.0-0
## [31] survival_3.2-11 zoo_1.8-9
## [33] glue_1.6.2 polyclip_1.10-0
## [35] gtable_0.3.0 zlibbioc_1.38.0
## [37] XVector_0.32.0 leiden_0.3.9
## [39] DelayedArray_0.18.0 future.apply_1.7.0
## [41] BiocGenerics_0.38.0 abind_1.4-5
## [43] scales_1.1.1 pheatmap_1.0.12
## [45] DBI_1.1.1 miniUI_0.1.1.1
## [47] Rcpp_1.0.7 viridisLite_0.4.0
## [49] xtable_1.8-4 reticulate_1.22
## [51] spatstat.core_2.1-2 bit_4.0.4
## [53] stats4_4.1.0 htmlwidgets_1.5.3
## [55] httr_1.4.2 RColorBrewer_1.1-2
## [57] ellipsis_0.3.2 ica_1.0-2
## [59] pkgconfig_2.0.3 farver_2.1.0
## [61] sass_0.4.0 uwot_0.1.10
## [63] deldir_1.0-6 utf8_1.2.1
## [65] tidyselect_1.2.0 labeling_0.4.2
## [67] rlang_1.0.6 reshape2_1.4.4
## [69] later_1.2.0 AnnotationDbi_1.54.0
## [71] munsell_0.5.0 tools_4.1.0
## [73] cachem_1.0.5 cli_3.4.1
## [75] generics_0.1.3 RSQLite_2.2.7
## [77] ggridges_0.5.3 org.Gg.eg.db_3.13.0
## [79] evaluate_0.20 fastmap_1.1.0
## [81] yaml_2.2.1 goftest_1.2-2
## [83] knitr_1.41 bit64_4.0.5
## [85] fitdistrplus_1.1-6 purrr_0.3.4
## [87] RANN_2.6.1 KEGGREST_1.32.0
## [89] pbapply_1.4-3 future_1.30.0
## [91] nlme_3.1-152 mime_0.10
## [93] compiler_4.1.0 rstudioapi_0.13
## [95] png_0.1-7 spatstat.utils_3.0-1
## [97] bslib_0.2.5.1 stringi_1.6.2
## [99] highr_0.9 lattice_0.20-44
## [101] Matrix_1.3-3 vctrs_0.5.1
## [103] pillar_1.8.1 lifecycle_1.0.3
## [105] spatstat.geom_3.0-3 lmtest_0.9-38
## [107] jquerylib_0.1.4 RcppAnnoy_0.0.19
## [109] bitops_1.0-7 data.table_1.14.0
## [111] cowplot_1.1.1 irlba_2.3.3
## [113] GenomicRanges_1.44.0 httpuv_1.6.1
## [115] R6_2.5.0 promises_1.2.0.1
## [117] KernSmooth_2.23-20 gridExtra_2.3
## [119] IRanges_2.26.0 parallelly_1.33.0
## [121] codetools_0.2-18 MASS_7.3-54
## [123] assertthat_0.2.1 SummarizedExperiment_1.22.0
## [125] withr_2.4.2 sctransform_0.3.3
## [127] GenomeInfoDbData_1.2.6 S4Vectors_0.30.0
## [129] mgcv_1.8-35 parallel_4.1.0
## [131] grid_4.1.0 rpart_4.1-15
## [133] tidyr_1.1.3 rmarkdown_2.17
## [135] MatrixGenerics_1.4.0 Rtsne_0.15
## [137] Biobase_2.52.0 shiny_1.6.0